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ABSTRACT 


3D  linear  elastic  finite  difference  calculations  were  used  to  investigate 
effects  of  near-source  scattering  on  long-period  surface  waves  radiated  by 
shallow  explosions  located  in  the  arc  of  a  subduction  zone.  These  simulations 
were  motivated  by  the  teleseismic  Rayleigh-wave  amplitude  patterns  of 
explosions  detonated  at  the  Amchitka  test  site.  Amplitude  patterns  of  these 
explosion-generated  20-40  second  surface  waves  show  evidence  of  both 
tectonic  release  and  source-receiver  path  effects. 

A  3D  velocity  model  was  constructed  based  on  the  Aleutian  subduction 
zone  structure  of  Boyd  and  Creager  (1991)  and  crustal  refraction  results 
reported  by  Lambert,  etai  (1970).  Two  3D  finite  difference  calculations  were 
performed  to  assess  the  sensitivity  of  the  results  to  the  source  location  and  fine 
details  of  the  velocity  model.  Each  finite  difference  grid  consisted  of  2  million 
nodes  and  the  calculations  were  designed  to  model  surface  waves  in  the  20  to 
40  second  period  range  at  distances  up  to  400  km  from  the  source. 

Displacements  on  the  free  surface  were  analyzed  to  show  the  effects  of 
the  velocity  structure  upon  the  3D  propagation  of  Rayleigh  waves  from  an 
explosion  source.  Azimuthaliy  dependent  Rayleigh  wave  amplitudes  are 
clearly  seen  as  well  as  Rayleigh-to-Love  wave  conversion  along  the  strike  of 
the  subduction  zone  structures.  Teleseismic  Rayleigh  wave  amplitude 
anomalies  of  up  to  a  factor  of  two  are  dearly  evident  from  the  near-source 
scattering.  Love  waves  1/3  to  1/4  the  amplitude  of  the  Rayleigh  waves  are 
radiated  in  narrow  azimuthal  ranges. 

A  method  for  the  continued  propagation  of  surface  waves  is  presented 
based  on  a  2D  Fresnel-Kirchoff  integral  for  surface  waves.  The  Fresnel-Kirchoff 
integral  accounts  for  diffraction  effects  assuming  uniform  propagation  outside 
the  finite  difference  grid.  This  hybrid  procedure  using  both  the  finite  difference 
simulation  and  the  Fresnel-Kirchoff  integral  predicts  50%  Rayleigh  wave 
amplitude  anomalies  at  teleseismic  distances  due  to  structure  within  400  km  of 
the  Amchitka  test  site. 


I.  INTRODUCTION 


Detailed  aspects  of  the  Rayleigh  wave  radiation  from  the  three  Amchitka 
explosions,  LONGSHOT,  MILROW,  and  CANNIKIN,  have  remained  a  problem 
since  the  early  1970's.  The  Rayleigh-wave  amplitudes  from  these  events  show 
evidence  for  both  tectonic  release  and  nonisotropic  source-receiver 
propagation  effects.  The  evidence  is  clear  that  the  large  events  were 
accompanied  by  tectonic  release.  However,  there  has  remained  questions  as 
to  how  much  of  the  apparent  radiation  pattern  from  the  three  Amchitka 
explosions  is  due  to  tectonic  release  and  how  much  is  due  to  source-receiver 
propagation  effects.  Since  Amchitka  Island  lies  above  an  active  subduction 
zone,  the  possibility  has  remained  that  the  complex  near-source  structure 
introduced  an  apparent  radiation  pattern  observed  at  teleseismic  distances.  In 
this  study,  we  model  the  long-period  surface-wave  propagation  in  the  near¬ 
source  subduction  zone  structure  around  Amchitka  Island  and  examine  the 
potential  nonisotropic  propagation  effects  on  Rayleigh  and  Love  waves  from 
isotropic  explosion  sources.  Three  dimensional  linear  elastic  finite  difference 
calculations  were  performed  to  model  the  complete  wave  propagation  effects 
for  periods  longer  than  20  seconds.  The  results  from  these  calculations  are 
compared  to  observed  Rayleigh  and  Love  waves  from  the  three  Amchitka 
explosions. 


2.  TECTONIC  RELEASE  VERSUS  PROPAGATION  EFFECTS 


In  order  to  reduce  systematic  errors  due  to  path  structure  and  put  the 
estimation  of  long-period  underground  explosion  source  characteristics  in  the 
form  of  a  a  physically  meaningful  quantity  such  as  the  seismic  moment  tensor, 
the  method  of  Rayleigh  wave  path  correction  was  introduced  (Stevens,  1986a). 
Rayleigh  and  Love  waves  from  the  test  site  are  used  to  derive  an  average 
structure  between  the  test  site  and  each  station  in  the  network.  The  Rayleigh 
and  Love  wave  amplitude  and  phase  (normal  or  reversed)  is  then  used  to 
estimate  a  moment  tensor  for  the  explosion  plus  tectonic  release  (Given  and 
Mailman,  1986;  Stevens,  1986a).  Unfortunately,  the  path  corrections  can  only 
partially  correct  for  the  effects  of  propagation  between  the  source  and  receivers. 
There  remain  propagation  effects  due  to  the  laterally  heterogeneous  structure  of 
the  Earth.  These  remaining  propagation  effects  result  in  individual  station 
corrections.  It  has  been  found  that  a  trade-off  exists  between  the  determination 
of  a  tectonic  release  pattern  for  events  at  the  test  site  and  the  average  station 
correction  pattern  for  the  network  of  stations.  Residual  propagation  effects  can 
be  mistaken  for  a  tectonic  release  pattern  and  vice  versa. 

Figure  1  shows  the  average  station  corrections  determined  from  the  three 
Amchitka  events  assuming  that  the  log  of  all  station  correction  factors  sum  to 
zero  (Stevens  and  McLaughlin,  1989).  Rather  than  showing  a  random 
azimuthal  pattern,  the  station  corrections  exhibit  a  strong  trend  from  large 
amplitudes  to  small  amplitudes  between  azimuths  of  40  and  70  degrees.  Rts  to 
the  Rayleigh  wave  amplitudes  for  CANNIKIN,  MILROW,  and  LONGSHOT  are 
shown  in  Figure  2  from  Stevens  and  McLaughlin  (1989)  with  F  factors  (ratio  of 
isotropic  to  double-couple  moment)  ranging  between  0.2  and  0.3.  Note  that  the 
inferred  tectonic  release  patterns  are  strongly  influenced  by  the  low  amplitudes 
to  the  east  and  large  amplitudes  to  the  north-east  as  seen  (inversely)  in  the 
station  corrections. 

This  azimuthal  radiation  pattern  has  been  studied  with  interest  since  the 
earliest  analysis  of  surface  waves  from  the  Amchitka  test  site.  Researchers 
have  compared  Rayleigh  waves  from  different  explosions  as  well  as  Rayleigh  to 
Love  wave  ratios,  and  studied  Rayleigh  waves  from  explosion  cavity  collapses 
in  order  to  estimate  the  effects  of  propagation  on  the  Rayleigh  wave  amplitudes 
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Figure  1 .  Path-corrected  long-period  Rayleigh-wave  (LR)  amplitudes  from 
Stevens  and  McLaughlin  (1989)  for  explosions  LONGSHOT, 
MILROW,  and  CANNIKIN  located  at  the  Amchitka  test  site.  Note 
the  trend  from  large  amplitudes  between  40  and  60  degrees  to 
small  amplitudes  between  60  and  90  degrees.  A  nonparametric 
smoother  has  been  applied  to  the  data  (dashed  line). 
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Amchitka  Radiation  Patterns 


Cannikin 

LoqMq  =  17.83  *  0.18 
F  =  0.21 
Strike  =  7* 


Miirow 

Log  Mq  =  17.14  *  0.21 
F  »  0.27 
Strike  =  337* 


Longshot 

Log  Mq  =  16.08  *  0.29 
F  =  0.30 
Strike  =  326* 


Figure  2.  Radiation  patterns  for  CANNIKIN  (top),  MILROW  (middle),  and 
LONGSHOT  (bottom)  inferred  from  path  corrected  Rayleigh  wave 
amplitudes.  The  trend  of  large  amplitudes  in  the  north-east  to 
small  amplitudes  in  the  east  result  in  radiation  patterns  with  a 
minimum  to  the  east-southeast  for  each  event.  The  ratios  of 
double-couple  to  isotropic  moment  range  from  0.2  to  0.3  for  the 
three  events. 
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used  to  estimate  the  surface  wave  magnitude.  (Ms)  and  moments  of  these 
explosions. 

Willis,  et  al.  (1972)  and  Toksoz  and  Kehrer  (1972)  both  studied  the 
Rayleigh  and  Lave  waves  radiated  by  CANNIKIN.  Toksoz  and  Kehrer 
concluded  from  an  analysis  of  Rayieigh/Love  wave  ratios  that  CANNIKIN  was 
accompanied  by  tectonic  release  (F  =  0.4)  with  an  isotropic  moment  of 
5  X  10^^  Nt-m.  Both  studies  noted  the  large  Rayleigh  wave  amplitudes  across 
the  Canadian  network  (azimuths  10  to  70  deg),  which  are  only  partly  explained 
by  a  tectonic  release  model.  In  addition,  Willis,  et  al.  examined  Rayleigh  wave 
amplitudes  from  the  MILROW  explosion  and  the  CANNIKIN  collapse.  The 
individual  MILROW  and  CANNIKIN  Ms  values  in  the  first  quadrant  (0-90  deg) 
from  Willis,  et  al.  are  shown  in  Rgure  3.  There  is  a  strong  gradient  in  Rayleigh 
wave  amplitude  of  nearly  0.5  magnitude  units  between  40  and  70  degrees.  The 
radiation  pattern  proposed  by  Toksoz  and  Kehrer  from  an  explosion  plus 
double-couple  Rayleigh-wave  radiation  pattern  is  shown  for  comparison  (dotted 
line).  Note  the  similarity  of  the  apparent  radiation  pattern  between  the  two 
events  suggesting  a  common  cause(s)  for  the  amplitude  patterns. 

In  Figure  4,  we  show  the  CANNIKIN  Rayleigh  wave  amplitudes  (from 
Willis,  et  al.)  normalized  to  Rayleigh  wave  amplitudes  from  the  CANNIKIN 
collapse  event.  It  is  often  assumed  that  a  coliapse  is  an  axisymmetric  source  so 
the  ratio  of  explosion  to  collapse  amplitudes  should  cancel  propagation  and 
instrumental  effects  that  otherwise  complicate  interpretation  of  Rayleigh  wave 
amplitudes  (Masse,  1971).  The  normalized  Rayleigh  waves  show  a  distinct 
azimuthal  pattern  that  is  partially  fit  by  the  theoreticai  radiation  pattern  of  an 
explosion  plus  a  double-couple  due  to  Toksoz  and  Kehrer. 

Also  comparing  explosions  and  collapses,  von  Seggern  (1973) 
concluded  from  a  study  of  MILROW,  MILROW  collapse,  and  LONGSHOT  Love 
waves  that  a  mechanism(s)  other  than  tectonic  release  contributed  to  Love 
waves  for  these  events.  He  used  matched  filters  derived  from  MILROWs  Love 
waves  to  isolate  the  Love  waves  excited  by  MILROW  collapse  and  the 
LONGSHOT  explosion.  Assuming  the  collapse  was  axisymmetnc,  it  should  not 
excite  Love  waves.  Furthermore,  the  ratio  of  Love  wave  amplitudes  between 
MILROW  and  LONGSHOT  are  consistent  with  the  yield  scaling  between  the  two 
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MILROW  &  CANNIKIN  Ms 


Azimuth  fDee) 


Figure  3.  Ms  station  magnitudes  for  CANNIKIN  (circles)  and  MILROW 
(triangles)  from  Willis,  et  al.  (1972).  The  radiation  pattern  from 
Toksoz  and  Kehrer  (1972)  is  superimposed  on  the  two  data  sets 
(dashed  line).  Note  the  trend  in  large  to  small  magnitudes 
between  40  and  90  degrees  azimuth. 
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Figure  4.  Normalized  station  amplitudes  for  CANNIKIN  from  Willis,  et  al. 

(1972).  The  radiation  pattern  from  Toksoz  and  Kehrer  (1972)  is 
superimposed  (dashed  line).  Note  the  trend  in  large  to  small 
magnitudes  between  40  and  90  degrees  azimuth  is  only  partially 
modeled  by  the  tectonic  release  pattern. 
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events,  suggesting  that  Love  wave  excitation  was  proportional  to  the  explosion 
part  of  the  source.  If  the  collapses  are  axisymmetric  sources,  then  von 
Seggern's  observations  are  consistent  with  the  hypothesis  that  the  MILROW 
and  LONGSHOT  Love  waves  were  produced  by  scattering  of  Rayleigh  waves  to 
Love  waves  along  the  propagation  path. 

Figure  5  shows  the  MILROW  explosion  and  collapse  Ms  magnitudes  from 
von  Seggern  and  Lambert  (1972).  The  amplitude  pattern  is  similar  to  Ms  data 
of  Willis,  et  al.  in  Figure  4  for  both  the  explosion  and  the  collapse.  Figure  6 
shows  the  MILROW  Ms  difference  (explosion  -  collapse)  at  common  stations  in 
the  first  quadrant.  No  statistically  significant  trend  is  present  in  the  Ms  difference 
between  40  and  90  degrees.  This  evidence  suggests  that  the  MILROW 
explosion  and  collapse  have  a  similar  radiation  pattern.  If  the  MILROW  collapse 
is  an  isotropic  source  then  it  appears  that  most  of  the  MILROW  explosion 
radiation  pattern  is  path  related  and  not  tectonic  release. 

The  results  of  von  Seggern’s  analysis  suggest  strong  near-source 
propagation  effects  that  affect  both  explosion  and  collapse  Rayleigh  and  Love 
excitation.  Since  the  subduction  zone  structure  near  Amchitka  Island  has  a 
reflection  symmetry  about  a  north-south  plane,  a  scattering  hypothesis  would 
predict  a  reflection  symmetry  in  the  Rayleigh  wave  radiation  pattern  about  a 
north-south  axis.  The  maximum  likelihood  WWSSN  Ms  station  corrections  from 
McLaughlin,  et  al.  (1986)  are  plotted  in  Figure  7  as  a  function  of  azimuth.  This 
study  utilized  a  large  number  of  WWSSN  stations  for  all  three  Amchitka  Island 
events  at  azimuths  not  used  in  some  other  studies.  The  station  effects  are 
plotted  from  0  to  360  deg  azimuth  in  the  lower  panel.  Note  the  dip  in  station 
effects  near  270  degrees  is  similar  to  the  dip  near  90  degrees.  In  the  upper 
panel,  stations  between  180  and  360  degrees  are  plotted  at  360  minus  the 
azimuth  as  would  be  suggested  by  a  reflection  symmetry  about  a  north-south 
line.  It  appears  from  Figure  7  that  the  WWSSN  Ms  station  effects  are  consistent 
with  a  north-south  reflection  symmetry  but  because  of  the  paucity  of  stations 
between  90  and  180  degrees  it  is  not  a  completely  convincing  test. 

It  is  interesting  that  based  on  data  from  collapse  events,  Willis,  et  al. 
(1972)  concluded  that  CANNIKIN  is  contaminated  by  tectonic  release  while  von 
Seggern  concluded  that  the  structural  effect  rather  than  tectonic  release  is  the 
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Figure  5.  MILROW  explosion  and  collapse  Ms  station  magnitudes  from  von 
Seggern  and  Lambert  (1972).  The  explosion  (triangles)  and 
collapse  (circles)  magnitudes  show  the  same  trend  between  40 
and  90  d^rees  azimuth.  Data  from  the  collapse  is  limited  to  the 
first  quadrant. 
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Figures.  MILROW  station  explosion  minus  collapse  Ms.  There  is  no 
significant  trend  in  the  data  between  40  and  90  degrees.  If  the 
collapse  was  an  axisymmetric  source  then  much  of  the  amplitude 
pattern  in  Figure  5  is  due  to  propagation  effects  and  not  tectonic 
release. 
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Figure  7.  Maximum  likelihood  station  effects  for  Ms  magnitudes  read  at 
WWSSN  stations  from  LONGSHOT.  MILROW,  and  CANNIKIN. 
Station  effects  are  plotted  from  0  to  360  degrees  (bottom)  and  from 
0  to  180  degrees  (top).  The  data  appears  to  be  consistent  with  a 
reflection  symmetry  about  a  north-south  axis  as  would  be 
suggested  by  the  symmetry  of  the  Aleutian  subduction  zone  in  the 
region  of  Amchitka.  A  nonparametric  smoother  has  been  applied 
to  the  data  (dashed  lines). 
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dominant  source  of  Love  waves.  Both  the  analysis  of  Willis,  et  al.  and  that  of 
von  Seggern  assumed  that  collapses  are  axisymmetric  sources,  yet  these 
authors  arrive  at  contradictory  conclusions.  Engdahl  (1972)  however,  presents 
double-couple  focal  mechanisms  for  CANNIKIN  collapse  events  that  suggest 
normal  faulting  on  45  degree  dipping  planes  striking  roughly  east-west.  These 
focal  mechanisms  are  based  on  locally  recorded  P-wave  first  motions.  It  is 
possible  that  these  double-couple  focal  mechanisms  do  not  reflect  the  actual 
long-period  nature  of  the  collapse;  however,  if  the  collapses  are  non-isotropic, 
then  the  simple  ratio  of  the  explosion-to-coilapse  does  not  represent  the  pure 
explosion  plus  double-couple  radiation  pattern.  Further  evidence  that  the  local 
stress  regime  is  influenced  by  the  collapse  is  the  observation  that 
microearthquake  activity  around  the  explosion  site  undergoes  an  abrupt 
reduction  immediately  following  the  collapse.  This  behavior  was  observed  for 
MILROW  and  CANNIKIN  (Engdahl,  1972)  and  is  common  for  events  at  NTS. 
Consequently,  although  the  explosion-to-collapse  amplitude  ratios  do  cancel 
propagation  and  instrumental  effects,  it  is  not  clear  that  they  reveal  a  pristine 
source  radiation  pattern  for  all  explosions.  The  resulting  pattern  may  reflect  the 
combined  source  radiation  pattern  of  both  the  explosion  and  collapse. 
Therefore  from  the  analyses  of  Willis,  et  al.  and  von  Seggern  it  is  dear  that  the 
explosions  are  contaminated  by  both  tectonic  release  and  path  propagation 
effects  but  by  an  undetermined  amount. 

In  addition  to  the  presence  of  Love  waves  and  the  strong  azimuthal 
radiation  patterns  from  Amchitka  explosions,  it  has  been  suggested  that  there 
was  an  apparent  bias  between  Rayleigh  wave  excitation  from  Amchitka 
explosions  and  NTS  explosions  of  similar  size,  von  Seggern  (1978)  examined 
Ms  for  megaton  NTS  events  compared  to  the  Ms  from  MILROW  and  CANNIKIN. 
He  concluded  that  Ms  from  the  Amchitka  events  was  low  relative  to  the  NTS 
events  by  as  much  as  0.5  magnitude  units.  Ms  values  determined  by  Marshall, 
et  al.  (1979)  and  moments  determined  by  Stevens  (1986b)  and  Stevens  and 
McLaughlin  (1989)  suggest  a  bias  closer  to  0.2  log-units. 

With  propagation  effects  superimposed  on  the  tectonic  release  radiation 
pattern,  it  is  possible  that  this  apparent  bias  is  due  to  the  nature  of  the  network 
sampling  afforded  by  the  WWSSN,  CSSN,  and  LRSM  seismic  networks  that 
were  used  to  estimate  Ms  and/or  moment  for  these  events.  The  networks  are 
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heavily  concentrated  in  North  America  and  Europe  where  azimuthally 
dependent  propagation  could  contribute  to  a  bias  in  the  estimation  of  the 
explosion  part  of  the  source.  The  most  noticeable  feature  of  this  azimuthal 
pattern  is  the  rapid  decrease  in  Rayleigh  wave  amplitude  with  change  in 
source-receiver  azimuth  between  40  and  70  degrees.  This  corresponds  to 
high  amplitudes  for  Canadian  stations  and  low  amplitudes  for  stations  in 
southern  North  America.  Coincidentally,  this  azimuthal  range  corresponds  to 
ray-paths  that  propagate  obliquely  across  the  Aleutian  subduction  zone  and  are 
obliquely  incident  upon  the  North  American  continental  margin.  The 
observation  that  average  WWSSN  Ms  station  effects  are  consistent  with  a 
reflection  about  a  north-south  plane  suggests  that  the  near-source  subduction 
zone  striking  east-west  has  introduced  an  apparent  radiation  pattern  into  the 
teleseismic  Rayleigh  waves. 

in  this  study  we  examine  the  hypothesis  that  near-source  3D  subduction 
structure  has  imposed  an  azimuthal  pattern  on  the  teleseismic  Rayleigh  wave 
amplitudes  from  the  Amchitka  test  site.  A  3D  elastic  finite  difference  code  is 
used  to  simulate  the  propagation  of  long-period  waves  in  a  model  of  the 
Aleutian  arc  near  Amchitka  Island. 


3.  3D  FINITE  DIFFERENCE  METHODS 


The  use  of  2D  finite  differences  has  become  increasingly  popular  to 
simulate  scalar  (acoustic)  and  elastodynamic  wavefields.  As  the  speed  and 
memory  capacity  of  computers  has  increased,  realistic  seismological  3D  finite 
difference  calculations  have  become  possible. 

The  principal  advantage  of  3D  elastodynamic  finite  differences  is  that  the 
calculations  implicitly  include  all  conversions  and  scattering.  The  principal 
disadvantage  is  that  a  limited  bandwidth  can  be  studied  and  the  computer  time 
is  proportional  to  the  fourth  power  of  the  size  of  the  model  and  bandwidth.  By 
limiting  the  bandwidth  ana  model  size,  long-period  near-source  scattering 
presents  a  tractable  problem  that  requires  only  a  few  compromises.  However, 
by  going  to  smaller  models,  the  bandwidth  of  interest  can  be  increased.  For 
example,  while  our  simulations  explore  0.05  Hz  waves  in  a  model  800  km 
across,  a  model  40  km  across  could  be  used  to  explore  1 .0  Hz  waves  with  the 
same  number  of  nodes.  Although  the  simulations  reported  in  this  paper  were 
performed  on  a  CRAY-2  supercomputer  using  vector  processing,  problems  of 
this  size  are  within  reach  of  many  modern  mini-supercomputers. 
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4.  THE  MODEL 


We  used  a  30  explicitly  time-stepping  linear  elastic  finite  difference 
method  to  calculate  the  response  to  an  explosion  source  of  a  model  for  the 
Aleutian  subduction  zone  in  the  region  of  Amchitka.  Our  objective  was  to 
accurately  model  elastic  waves  with  periods  20  seconds  or  longer  over  a  region 
extending  400  km  from  Amchitka  island  in  all  directions.  To  achieve  this 
resolution  and  spatial  extent  efficiently,  we  generated  a  finite  difference  grid 
with  an  inner,  finely  gridded  region  with  6  km  spacing.  Outside  this  high- 
reiolution  inner  grid,  the  spacing  was  slowly  expanded  at  a  constant  rate  of 
increase.  In  the  region  of  expanding  grid,  artificial  viscosity  was  increased  to 
damp  reflections  from  the  expanding  grid  and  the  outer  boundaries  of  the  3D 
grid.  In  this  way,  reflections  from  the  outer  edges  are  eliminated  and  do  not 
contaminate  the  inner  portion  of  the  grid,  in  addition,  we  took  advantage  of  a 
natural  symmetry  in  the  problem,  assuming  reflection  symmetry  about  a  north- 
south  vertical  plane  passing  through  the  source.  The  inner  portion  of  the  grid 
was  70x70x140  grid  points  (420  by  420  by  840  km).  The  time  step  was  0.25 
sec.  Figure  8  illustrates  this  use  of  an  expanding,  attenuating  grid  outside  the 
inner  grid.  The  total  grid  used  in  the  computation  was  100  by  100  by  200  or 
2  million  grid  points.  A  840  by  840  km  box  is  superposed  on  top  of  a  map  of 
the  Amchitka  region  in  Rgure  9.  The  P-wave  velocity,  S-wave  velocity,  density, 
and  viscosity  were  specified  at  each  grid  point. 

The  velocity  structure  for  the  model  was  derived  from  two  sources. 
Shallow  structure  for  the  area  from  marine  seismic  profiles  was  reviewed  by 
Lambert,  etal.  (1970).  They  present  a  model  with  maximum  crustal  thickness 
between  the  island  and  the  trench  (about  40  km),  and  thinning  crust  behind  the 
arc  to  merge  with  a  normal  oceanic  crust.  We  chose  the  oceanic  PREM  model 
north  and  south  of  the  subduction  zone  as  the  limiting  far-fieid  velocity  model. 
The  mantle  velocity  from  the  spherically  symmetric  PREM  model  (Dziewonski 
and  Anderson,  1983)  was  perturbed  in  accordance  with  a  thermal  model  for  the 
subducting  lithosphere  from  Boyd  and  Creager  (1991).  Boyd  and  Creager  used 
local  seismicity,  P-wave  residuals  sphere  analysis  and  travel  times  from  the 
CANNIKIN  explosion  to  constrain  the  geometry  of  the  subducting  Pacific  plate 
under  the  Aleutian  Islands.  Their  analysis  was  guided  by  a  thermal  model  for 
the  subducting  slab.  The  study  of  Boyd  and  Creager  follows  a  long  line  of 
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3-D  Finite  Difference  Grid  Expanding  and  Attenuating  Grid 


Figure  8.  The  30  finite  difference  grid  was  oriented  with  the  v  axis  aligned  north,  u  axis  east  and  z  axis  down.  A 
refit  . tion  symmetry  boundary  condition  was  used  on  the  u  =  0  plane  taking  advantage  of  the  north-south 
reflection  symmetry  of  the  subduction  zone.  The  inner  grid  had  6  km  spacing  and  was  1 40x70x70.  The 
outer  portion  of  the  grid  had  a  progressively  expanding  spacing.  The  entire  grid  was  200x100x100. 
Numerical  viscosity  was  used  to  damp  reflections  from  the  expanding  grid.  A  0.25  second  time  step  was 
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Figure  9. 


Map  of  region  around  Amchitka  Island  test  site  showing  the 
approximate  location  of  the  trench  axis  and  an  840  by  840  km  box 
representing  the  finite  difference  gn'd  used  in  the  simulations. 
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studies  that  have  examined  teleseismic  P-wave  travel  time  anomalies  from 
Amchitka  Island  explosions  (for  example  see  Chiburis  and  Ahner,  1969;  Jacob, 
1972;  Page,  1972;  Davies  and  Julian,  1972;  Sleep,  1973;  Cormier,  1990). 
P-wave  and  S-wave  velocity  contours  are  shown  for  the  finite  difference  model 
in  Figure  10. 

The  oceanic  water  layer  was  not  included  in  the  3D  calculations.  The 
effect  of  the  island/ocean  interface  was  instead  modeled  in  a  separate  series  of 
2D  axisymmetric  finite  difference  simulations  (Stevens,  etal.,  1990).  The  results 
of  this  study  show  that  in  cases  where  an  island  rises  steeply  out  of  the  ocean, 
the  Rayleigh  waves  excited  by  an  explosion  can  be  sharply  reduced.  However, 
for  the  specific  case  of  Amchitka  Island,  the  bathymetry  is  too  gentle  to  cause 
any  significant  amplitude  reduction;  the  effect  of  the  ocean  and  island/ocean 
interface  on  20-50  second  Rayleigh  waves  is  negligible. 

Two  3D  linear  elastic  finite  difference  calculations  were  performed.  In  the 
first  calculation,  curvature  of  the  Aleutian  arc  was  ignored  and  the  source  was  at 
location  #1  indicated  in  Figure  10.  In  the  second  calculation,  curvature  of  the 
subduction  zone  structure  was  included  and  the  source  location  is  indicated  by 
#2  in  Figure  10.  Locations  #1  and  #2  bracket  the  location  of  Amchitka  Island  in 
the  subduction  zone  structure,  with  the  first  closest  to  the  approximate  position 
of  the  island.  A  comparison  of  the  two  calculations  shows  similarities  in  the 
azimuthal  effects,  but  it  is  clear  that  the  final  results  are  sensitive  to  the  location 
of  the  source  in  the  model  and  some  model  details.  Some  aspects  of  the 
velocity  model  are  poorly  controlled,  such  as  the  detailed  thickness  of  the  crust 
and  the  low  velocity  accretion  zone  north  of  the  trench.  For  this  reason,  results 
are  best  considered  as  indicative  of  the  magnitude  and  qualitative  character  of 
the  probable  effects  of  the  subduction  zone  upon  wave  propagation. 


Figure  1 0.  Contours  of  the  S-wave  velocity  field  along  the  north-south  vertical 
plane.  The  model  is  based  on  Boyd  and  Creager  (1991)  and 
Lambert,  et  al  (1 970).  Two  simulations  were  conducted,  with  an 
explosive  source  at  #1  and  #2  indicated  in  the  figure. 
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5.  ANALYSIS  AND  DISCUSSION 


Several  analysis  procedures  were  applied  to  the  displacements  on  the 
surface  of  the  grid.  These  included  animated  snap-shots  of  the  displacement 
field  at  selected  times,  as  well  as  procedures  adapted  from  traditional  surface 
wave  data  analysis,  narrow  band  filter  amplitude  estimation,  group  and  phase 
velocity  estimation,  and  frequency-wavenumber  analysis. 

Animated  snap-shots  of  the  displacement  field  were  resolved  into 
vertical,  radial,  and  transverse  components  in  order  to  visualize  the  wave 
propagation.  A  color  video  of  the  evolving  wavefield  was  made.  The  video 
movie  was  most  helpful  in  visualizing  the  wave  propagation  and  identifying 
areas  for  more  detailed  quantitative  study. 

Several  monochrome  (contour  plot)  snap-shots  are  shown  in  Figures 
1 1 A  through  1 1 D.  The  asymmetries  in  group  velocity  and  amplitude  are  clearly 
evident  from  these  figures.  The  structure  to  the  east  of  the  source  along  the 
strike  of  the  structure  is  slower  in  the  20-30  second  period  range  because  of  the 
low  velocity  sediments  and  thick  crust  along  the  trench  axis.  Consequently, 
there  exists  a  narrow  low-velocity  zone  that  forms  a  channel  for  surface  waves. 
Refraction  around  this  channel  causes  focusing  of  surface  wave  amplitudes 
along  the  axis  of  the  channel  and  defocusing  in  regions  immediately  off-axis 
(north  and  south)  of  the  channel. 

The  snap-shots  (Figure  1 1 )  of  the  transverse  component  of  motion  show 
a  growing  transverse  wave  that  develops  along  the  north  and  south  margins  of 
the  low  velocity  channel.  The  transverse  amplitudes  are  about  1/4  the 
amplitude  of  the  radial  Rayleigh  wave  and  appear  to  be  associated  with  the 
passage  of  the  Rayieigh  wave  along  the  edges  of  the  low-velocity  channel. 
Frequency-wavenumber  spectra  were  used  to  confirm  the  presence  of  Love 
wave  type  motion. 

Arrays  of  grid  points  were  selected  at  locations  on  the  surface  of  the  grid 
and  frequency-wavenumber  (FK)  spectra  were  estimated  from  the  three 
components  of  motion.  Figures  12  and  13  show  FK  spectra  from  synthetic 
arrays  at  azimuths  of  68  and  160  degrees  respectively.  Array  #1  at  an  azimuth 
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Figure  11  A.  Contours  of  displacement  amplitude  of  the  vertical,  radial,  and 
transverse  components  of  motion  on  the  free-surface  of  the  finite 
difference  grid  at  time  t »  70  seconds.  The  motion  is  dominated  by 
the  fundamental  Rayleigh  wave  on  the  vertical  and  radial 
components.  Note  the  transverse  motion  beginning  to  form  behind 
the  Rayleigh  wave  southeast  and  northeast  of  the  source. 
Simulation  30#1. 


21 


North  (km) 


Vertical.  t»100  sec.  Radial.  t«100  sec.  Transverse.  t«100  sec. 


last  (km)  Kast  (km)  East  (km) 


Figure  11 B.  Contours  of  displacement  amplitude  of  the  vertical,  radial,  and 
transverse  components  of  motion  on  the  free-surface  of  the  finite 
difference  grid  at  time  t  >  100  seconds.  Asymmetry  in  the  Rayleigh 
wave  amplitude  and  arrival  time  can  be  seen  in  the  vertical  and 
radial  components.  Note  the  transverse  Love  wave  motion  behind 
the  Rayleigh  wave  southeast  and  northeast  of  the  source. 
Simulation  3D#1. 
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Figure  11C.  Contours  of  displacement  amplitude  of  the  vertical,  radial,  and 
transverse  components  of  motion  on  the  free-surface  of  the  finite 
difference  grid  at  time  t  s  130  seconds.  Rayleigh  wave  motion  is 
both  stronger  and  delayed  to  the  east  along  the  axis  of  the 
subduction  zone.  Love  wave  motion  southeast  and  northeast  of 
the  source  is  beginning  to  overtake  the  Rayleigh  wave  motion. 
Simulation  3D#1. 
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Figure  1 1 D.  Contours  of  displacement  amplitude  of  the  vertical,  radial,  and 
transverse  components  of  motion  on  the  free-surface  of  the  finite 
difference  grid  at  time  t  *  1 50  seconds.  Rayleigh  wave  motion  is 
both  stronger  and  delayed  to  the  east  along  the  axis  of  the 
subduction  zone.  Love  wave  motion  southeast  and  northeast  of 
the  source  is  overtaking  the  Rayleigh  wave  motion.  Simulation 
3D#1. 
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Figure  1 2.  Wavenumber  spectra  at  0.047  Hz  from  a  synthetic  array  of  stations 

located  300  km  from  the  source  at  an  azimuth  of  68  degrees.  The 
vertical  and  east  components  show  a  prominent  arrival  with  a 
phase  velocity  of  about  3.3  km/sec  propagating  at  an  azimuth  of 
80  degrees.  The  north  component  shows  an  arrival  propagating  at 
4.0  km/sec  and  1 07  degrees.  The  vertical  and  east  component  see 
a  12  degree  off  azimuth  Rayleigh  wave  while  the  north  component 
records  a  40  degrees  off  azimuth  Love  wave. 
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Figure  1 3.  Wavenumber  spectra  at  0.047  Hz  from  a  synthetic  array  of  stations 
located  300  km  from  the  source  at  an  azimuth  of  160  degrees.  All 
three  components  show  a  prominent  arrival  with  a  phase  velocity 
of  3.44  km/sec  propagating  at  an  azimuth  of  158  degrees.  The 
Rayleigh  wave  is  on  azimuth  from  the  source. 
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of  68  degrees  shows  a  prominent  arrival  propagating  with  a  phase  velocity  of 
3.2-3.4  km/sec  at  an  azimuth  of  80  degrees  on  the  vertical  and  east 
components.  The  north  component  shows  a  prominent  arrival  propagating  with 
a  phase  velocity  of  4.0  km/sec  at  an  azimuth  of  107  degrees.  The  FK  spectra 
demonstrate  that  the  transverse  component  waves  (north)  propagated  at  Love 
wave  phase  velocities  but  in  a  direction  40  degrees  from  the  source  to  receiver 
azimuth.  All  three  components  of  FK  spectra  demonstrate  that  both  the 
Rayleigh  and  Love  wave  fields  are  refracted  along  the  low-velocity  channel 
associated  with  the  thick  crust  and  low  velocity  accretionary  wedge.  Array  #2 
located  at  an  azimuth  of  160  degrees  from  the  source  shows  ail  three 
components  dominated  by  a  Rayleigh  wave  propagating  at  a  phase  velocity  of 
3.44  km/s  and  an  azimuth  of  158  degrees.  No  Love  waves  are  seen  in  the 
direction  of  array  #2.  and  the  Rayleigh  wave  is  on  azimuth. 

In  order  to  quantify  the  amplitude  and  velocity  anomalies  associated  with 
the  30  wave  propagation,  conventional  narrow  bandpass  filtering  was  used  to 
examine  the  seismograms  available  across  the  free-surface  of  the  grid.  Figure 
14A  shows  the  vertical  component  traces  for  a  distance  of  270  km  from  the 
source.  The  traces  are  shown  from  0  to  180  degrees  azimuth.  Note  the  larger 
amplitudes  and  delayed  waveforms  to  the  east.  Similar  effects  are  seen  on  the 
radial  components,  shown  in  Figure  14B  at  the  same  scale  as  the  vertical 
components.  The  transverse  components  of  motion  are  displayed  in  Figure 
14C  at  a  larger  scale  (5X)  than  the  vertical/radial  traces.  The  transverse 
components  are  largest  in  the  60-80  and  100-130  degree  azimuth  ranges. 

Figure  15  shows  the  peak  bandpass  filtered  Rayleigh-wave  amplitudes 
as  a  function  of  azimuth  for  two  frequencies  and  both  simulations.  Note  that  the 
large  amplitudes  tend  to  occur  at  the  easterly  azimuths  but  that  the  pattern  is 
quite  different  for  the  two  source  locations  (3D#1  and  3D#2).  Figure  16  shows 
the  peak  Rayleigh-wave  amplitudes  as  a  function  of  azimuth  for  a  fixed  distance 
of  350  km  compared  to  the  arrival  time.  Note  that  the  larger  amplitudes  tend  to 
arrive  late.  Higher  amplitudes  tend  to  occur  in  slower  directions.  This  analysis 
supports  the  contention  that  the  low  velocity  region  along  the  subduction  zone 
focuses  the  surface  wave  energy.  Because  the  sources  are  offset  from  the 
center  of  the  wedge/trench  axis  the  patterns  are  not  symmetrical  about  the  east- 
west  azimuth. 
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Figure  14A.  Vertical  component  of  motion  displayed  for  a  constant  distance 
from  the  source,  270  km,  as  a  function  of  time  and  azimuth.  Note 
that  the  traces  are  delayed  and  larger  amplitude  to  the  east. 


Azimuth  (Deg) 


Azimuth  (Deg) 


Azimulh  (Dog) 


Relative  Amplitude  Relative  Amplitude 


I 


UR  Peak  Amplitude,  A  =  300km,  3D#1 

1.4 
1.2 
1.0 
0.8 
0.6 

0  40  80  120  160 


UR  Peak  Amplitude,  A  =  300km,  3D#2 

1.4 
1.2 
1.0 
0.8 
0.6 

0  40  80  120  160 

>  Azimuth  (deg) 


Figure  1 5.  Peak  bandpass  filtered  Rayleigh  wave  (LR)  amplitude  (periods  of 
20  and  30  sec)  as  a  function  of  azimuth  at  a  distance  of  300  km 
from  the  source.  Simulations  3D#1  and  3D#2  show  similar 
patterns  but  significant  differences  in  details  of  the  patterns.  Both 
simulations  show  factors  of  2  in  azimuthal  variation  at  300  km  from 
the  source. 
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Rgure  1 6.  Peak  20  and  30  second  Rayleigh  wave  amplitude  at  350  km  from 

the  source  (bottom)  and  the  arrival  time  (top)  for  simulation  3D#1 . 
Note  the  correlation  between  the  arrival  time  and  the  amplitude. 
Amplitudes  tend  to  be  higher  in  slower  directions.  Focusing 
occurs  along  the  trench  axis  70-100  degrees. 
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In  order  to  show  the  oetails  of  the  surface  wave  amplitude  anomalies 
across  the  surface  of  the  grid,  each  vertical  seismogram  on  the  grid  surface  was 
narrow  band  filtered  and  the  Rayleigh  wave  amplitude  was  estimated. 
Examples  of  these  amplitudes  are  contoured  in  Figures  17A  and  18A.  The 
Amplitude  contours  show  the  effects  of  both  normal  geometric  attenuation  and 
azimuthal  variation  from  the  3D  propagation.  Figures  17B  and  18B  show  the 
amplitude  contours  corrected  for  1/Vr  geometrical  spreading.  Due  to  near>field 
terms,  the  corrected  amplitudes  would  not  be  constant  even  in  the  absence  of 
lateral  heterogeneity,  but  the  1/Vr  correction  does  serve  to  bring  out  the 
azimuthal  anomalies  and  the  fact  that  the  anomalies  tend  to  grow  with  distance 
as  the  focusing/defocusing  lateral  refraction  progressively  disturbs  the  wave 
field. 


To  compute  Rayleigh  amplitudes  at  teleseismic  distances,  we  developed 
a  Fresnel-Kirchoff  integral  procedure  to  propagate  surface  waves  to  teleseismic 
distances  and  account  for  diffraction  of  the  wavefront.  A  derivation  of  the 
application  of  Fresnel-Kirchoff  diffraction  theory  to  surface  waves  is  contained  in 
Appendix  I.  This  method  assumes  that  the  surface  waves  propagate  from  the 
edge  of  the  3D  grid  into  a  laterally  homogeneous  earth  structure  to  receivers  far 
away  without  additional  complication.  This  assumption  allows  us  to  calculate 
the  effects  of  far-field  diffraction  from  the  near-source  scattering  modeled  by  the 
3D  finite  difference  grid. 

To  apply  the  Fresnel-Kirchoff  integral  we  save  the  synthetic  free-surface 
displacement  seismograms,  u(  j;',  t),  on  a  circle  of  constant  radius,  R,  from  the 
source  and  apply  the  frequency-domain  integral, 

k  ^7  fi'kfczd 

u(L  =  ^  J  u(r',<a)^-— (cos(0')  +  l)Rd0'  (9) 

where  0'=  0  -  0o  is  the  angle  between  the  far-field  receiver  and  the  seismogram 
Ircation. 
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Figure  17A.  Contours  of  peak  LR  amplitude  at  0.05  Hz  from  simulation  3D#1. 
Note  the  elongation  of  the  contours  in  the  east  west  direction. 
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Figure  17B.  Contours  of  peak  LR  amplitude  at  0.05  Hz  from  3D  simulation 

3D#1,  corrected  for  1/Vr  geometrical  spreading.  Note  the 
elongation  of  the  contours  in  the  east  west  direction. 
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Figure  18A.  Contours  of  peak  LR  amplitude  at  0.05  Hz  from  simulation  3D#2. 
Note  the  elongation  of  the  contours  in  the  east  west  direction. 


LR  Amplitudes,  3D#2,  T=20  sec 


<-WE->  (km) 
Amplitudes  *  Sqrt(R) 


Figure  186.  Contours  of  peak  LR  amplitude  at  0.05  Hz  from  simulation  3D#2, 

corrected  for  1/Vr  geometrical  spreading.  Note  the  elongation  of 
the  contours  in  the  east  west  direction. 
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The  amplitude  and  phase  of  the  free-surface  seismograms  at  fixed 
distances  from  the  source  are  plotted  in  Figures  19A  and  19B  for  the  two 
simulations.  This  information  is  utilized  in  the  numerical  integral  to  compute  the 
far-field  teleseismic  surface  wave  amplitudes.  The  phase  (delay-advance)  of 
the  wavefield  as  a  function  of  azimuth  at  the  integration  distance,  R,  is  just  as 
important  as  the  amplitude  dependence  in  the  final  amplitude  at  the  far-field 
receiver. 

Note  the  negative  correlation  between  amplitude  and  phase  for 
simulation  3D#1  in  Figure  19A.  This  suggests  that  focusing/defocusing  is 
largely  responsible  for  the  amplitude  pattern  at  this  distance  from  the  source. 
Note  also  that  for  both  simulations  the  amplitude  patterns  are  similar  at  different 
frequencies  while  the  detailed  patterns  are  frequency  dependent.  Figures  20A 
and  20B  plot  the  far<field  relative  amplitude  as  a  function  of  amplitude  and 
frequency.  Although  the  patterns  are  frequency  dependent,  the  different 
frequencies  share  some  common  patterns.  The  average  relative  pattern  from 
32  to  20  seconds  period  is  shown  in  the  upper  frame  of  Rgures  20  A  and  B. 

The  dominant  pattern  of  the  observed  amplitudes  was  high  amplitudes  in 
the  northeast  and  low  amplitudes  to  the  east  (Rgures  1-  7).  Neither  simulation, 
3D#1  or  3D#2,  reproduces  the  exact  observed  pattern,  but  both  simulations 
show  focusing/defocusing  patterns  that  have  characteristics  similar  to  the 
observed  amplitude/magnitude  patterns.  The  principal  characteristic  that  both 
data  and  simulation  exhibit  is  azimuth  ranges  of  low  amplitude  adjacent  to 
azimuth  ranges  of  high  amplitude.  Both  the  data  and  the  simulated  amplitude 
patterns  show  azimuthally  dependent  variations  that  can  exceed  factors  of  2 
peak-to-peak  over  a  broad  frequency  range.  The  amplitude  patterns  can  be 
very  narrow  (10  or  20  degrees),  uncharacteristic  of  the  broad  amplitude  patterns 
from  tectonic  release  which  have  amplitude  perturbations  proportional  to 
sin(2(0  -  0’)). 

Although  the  hybrid  finite-difference  Fresnel-Kirchoff  synthetic  amplitude 
patterns  do  not  reproduce  the  observed  amplitude  patterns  it  is  clear  that  they 
do  predict  variations  comparable  to  those  observed.  Furthermore,  the  two 
simulations,  3D#1  and  3D#2,  illustrate  that  changes  in  the  location  of  the  source 
comparable  to  a  wavelength  (32  seconds)  alter  the  amplitude  pattern 
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Phase  vs  Azimuth  at  300  Km,  3D#1 
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Figure  19A.  The  relative  amplitude  (txjttom)  and  phase  (top)  for  the  Rayleigh 
wave  at  a  distance  of  300  km  from  the  source  in  3D  simulation  #1 . 
Note  the  negative  correlation  between  the  large  amplitude  and 
phase  (delayed)  for  azimuths  to  the  east  or  west  (90  degrees). 
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Figure  19B.  The  relative  amplitude  (bottom)  and  phase  (top)  for  the  Rayleigh 
wave  at  a  distance  of  370  km  from  the  source  in  30  simulation  #2. 
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Amplitude  vs  Azimuth  at  10000  Km,  3D#1 


Figure  20A.  Predicted  teleseismic  amplitudes  from  30  simulation  #1  as  a 
function  of  azimuth  using  the  Fresnel-Kirchoff  integral  to  account 
for  diffraction  effects  in  the  far-fleld.  The  average  relative  amplitude 
from  32  to  20  seconds  period  is  at  the  top. 
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Amplitude  vs  Azimuth  at  10000  Km,  3D#2 
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Figure  206.  Predicted  teleseismic  amplitudes  from  3D  simulation  #2  as  a 
function  of  azimuth  using  the  Fresnei-Kirchoff  integral  to  account 
for  diffraction  effects  in  the  ^r-fieid.  The  average  relative  amplitude 
from  32  to  20  seconds  period  is  at  the  top. 
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significantly.  Similarly  the  results  appear  to  be  sensitive  to  details  of  the  model 
and  in  particular  the  thickness  of  the  crust  in  the  region  around  the  source.  The 
local  shallow  model  was  poorly  constrained. 
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6.  CONCLUSIONS 


We  present  results  from  a  large  3D  linear  elastic  finite  difference 
calculation  designed  to  simulate  20-40  second  surface  waves.  Furthermore,  we 
propose  a  procedure  based  on  the  Fresnel-Kirchoff  integral  for  surface  waves  to 
continue  the  propagation  from  a  3D  finite  differences  grid  to  teleseismic 
distances.  The  combination  of  these  two  techniques  allow  the  seismologist  to 
explore  3D  near-source  scattering  effects  on  surface  waves  and  project  the 
results  to  teleseismic  distances.  3D  finite  differences  is  a  complete  method  that 
includes  conversions,  while  the  Fresnel-Kirchoff  integral  assumes  a  layered 
earth  structure  outside  the  3D  grid  and  accounts  for  far-field  diffraction  effects. 

A  3D  model  for  the  velocity  structure  of  the  Aleutian  arc  has  been  used  to 
predict  teleseismic  surface  wave  amplitudes  from  a  shallow  explosive  source 
located  on  Amchitka  island.  The  predicted  Rayleigh  wave  amplitudes  show  as 
much  as  a  factor  of  two  variation  with  azimuth.  The  amount  of  azimuthal 
variation  is  consistent  with  observations,  although  the  observed  pattern  of  the 
variation  with  azimuth  is  not  reproduced  exactly.  Given  that  the  shallow  3D 
velocity  structure  of  the  arc  is  not  well  known,  it  is  encouraging  that  the  simple 
3D  structure  predicts  an  anomaly  of  about  the  right  magnitude  and  similar 
character  to  that  observed.  Furthermore,  the  3D  model  was  restricted  to  model 
heterogeneity  within  400  km  of  the  Amchitka  test  site,  and  surface  waves  from 
Amchitka  must  propagate  further  along  laterally  varying  structure  to  north 
American  stations.  It  appears  that  much  of  the  observed  surface  wave 
amplitude  patterns  from  Amchitka  island  can  be  explained  by  refraction  effects 
near  the  source  region.  The  3D  modeling  also  predicts  teleseismic  Love  waves 
from  Rayleigh  to  Love  scattering.  Love  waves  were  indeed  observed  by  von 
Seggern  (1973)  from  the  MILROW  collapse.  These  Love  waves  would  be 
observed  over  limited  azimuths,  and  with  amplitudes  approximately  1/3  to  1/4 
the  Rayleigh  wave  amplitude. 

The  dominant  mechanism  for  the  Rayleigh  wave  amplitude  variation  in 
the  model  appears  to  be  lateral  refraction  due  to  a  low  velocity  channel.  This 
low  velocity  channel  in  the  20  to  30  second  bandwidth  is  caused  by  a  thick  crust 
and  low  velocity  sediments  along  the  arc.  The  refraction  is  frequency 
dependent,  and,  because  of  the  high  velocity  slab  beneath  the  trench  axis,  this 
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region  is  a  high  velocity  zone  at  periods  longer  than  35  seconds.  Sources  in 
and  near  the  low  velocity  channel  will  appear  to  have  a  far-field  radiation 
pattern  at  teleseismic  distances  because  of  the  lateral  refraction  of  energy  along 
the  channel.  Similar  focusing/defocusing  effects  have  been  modeled  using  ray 
tracing  (von  Seggern,  atal.,  1975),  Gaussian  beams  (Zeng,  etal.,  1989)  and  a 
membrane  vibration  model  of  surface  waves  (Tanimoto,  1990).  Ray  tracing, 
Gaussian  beams,  and  membrane  models  for  surface  wave  refraction  do  not 
include  mode  conversions  observed  in  the  3D  elastic  modeling.  Rayleigh  to 
Love  wave  scattering  was  strongest  along  the  edges  of  the  low  velocity  channel 
where  lateral  gradients  are  maximum  and  it  is  reasonable  to  assume  that  the 
Rayleigh  to  Love  and  Love  to  Rayleigh  scattering  phenomenon  is  common  in 
regions  of  such  strong  lateral  contrasts. 

Assuming  that  near-source  focusing/defocusing  is  a  significant 
mechanism  for  the  azimuthally  dependent  amplitude  variations,  we  can  begin  to 
estimate  the  possible  bias  this  may  introduce  into  Ms  or  explosion  moment 
determinations.  The  rms  log-amplitude  azimuthal  variations  predicted  by  30 
simulations  #1  and  #2  (top  of  Figures  20  A  and  B)  are  both  0.1  magnitude  units. 
By  contrast  the  rms  station  corrections  of  Stevens  and  McLaughlin  (1989) 
(Figure  1)  and  McLaughlin,  etal.  (1986)  (Figure  7)  are  about  0.2  magnitude 
units.  The  rms  smoothed  azimuthal  variations  (dashed  lines  in  Figures  1  and  7) 
are  about  0.12  to  0.14.  Smoothing  removes  some  of  the  scatter  due  to 
additional  propagation  and  station  effects  far  from  the  near-source  region.  If  we 
assume  that  the  smooth  amplitude  versus  azimuth  effects  are  caused  by  near¬ 
source  focusing/defocusing  then  we  can  ask  what  possible  bias  may  have  been 
introduced  into  the  estimated  surface  wave  amplitude. 

There  are  two  possible  sources  of  estimation  bias.  The  first  source 
comes  from  the  non-uniform  azimuthal  sampling  of  data.  Given  the  gaps  in 
azimuth  for  which  there  is  no  or  little  data  this  source  of  bias  can  only  be  partly 
addressed  since  we  were  not  able  to  reproduce  the  observed  amplitude 
pattern.  However,  if  we  uniformly  re-sample  the  smoothed  amplitude  patterns  in 
Figures  1  and  7  there  appears  to  be  no  bias  due  to  the  non-uniform  sampling. 

The  second  source  of  bias  comes  from  log-amplitude  averaging  of  the 
focusing/defocusing  pattern  as  apposed  to  rms-amplitude  averaging.  By 
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definition  the  log-average  of  the  station  effects  are  zero  as  determined  by 
Stevens  and  McLaughlin  (1989)  or  McLaughlin,  etal.  (1986)  in  Figures  1  and  7 
respectively.  However,  the  use  of  log-averaging  inherent  in  the  statistical 
models  that  they  used  may  lead  to  a  bias  in  the  estimated  total  energy  in  the 
presence  of  focusing/defocusing.  if  the  smoothed  amplitude  variation  is  due  to 
focusing/defocusing  then  total  radiated  surface  wave  energy  integrated  over 
azimuth  should  be  invariant;  the  rms  amplitude  sampled  uniformly  over  all 
azimuths  should  be  unbiased  (McLaughlin,  1986).  Focusing/defocusing  only 
redistributes  energy  from  one  azimuthal  range  to  another  and  does  not  change 
the  total  energy  radiated  by  the  source.  Assuming  that  the  smoothed  amplitude 
function  is  from  near-source  focusing/defocusing  we  can  estimate  the  bias  by 
resampling  the  smoothed  amplitude  functions  and  compute  the  rms  amplitude. 
For  the  smoothed  amplitude  function  in  Figure  1,  the  uniformly  sampled  rms 
amplitude  is  1.64  (0.21  magnitude  units)  larger  than  the  baseline.  Similarly,  the 
estimated  bias  from  Figure  7  is  1.47  or  0.17  magnitude  units.  Consequently,  the 
log-average  of  surface  wave  amplitude  as  sampled  by  the  WWSSN  network  is 
about  0.17  to  0.21  log-units  smaller  than  it  would  be  if  the  assumed  near-source 
focusing/defocusing  (dashed  lines  in  Figures  1  and  7)  were  absent. 
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APPENDIX  I. 


USE  OF  THE  SURFACE-WAVE  FRESNEL-KIRCHOFF 
INTEGRAL  WITH  FINITE  DIFFERENCE  CALCULATIONS 
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With  a  few  justifiable  assumptions  we  can  compute  far-field  Rayleigh 
waves  from  the  monitored  displacements  on  a  closed  loop  located  on  the 
surface  of  the  3D  finite  difference  grid.  We  use  the  Fresnel-Kirchoff  integral  to 
compute  the  teleseismic  response  of  the  surface  waves  assuming  that  the 
region  outside  of  the  3D  grid  is  uniform  or  that  the  propagator  from  the  outer 
edge  of  the  3D  grid  can  be  written  simply  as  a  phase  factor  and  a  geometrical 
spreading  term. 

For  simplicity  assume  that  the  far-field  propagator  for  the  surface  wave 
outside  of  the  3D  finite  difference  grid  is 


P(|r-r1)= 


(1) 


and  that  the  surface  wave  amplitude  incident  of  frequency,  co,  incident  upon  a 
closed  loop  at  the  free  surface.  S,  is  given  by  u(r',  cb)  =  A(r‘,  o))e“l>(“),  then  the 
Fresnel-Kirchoff  integral  for  the  surface  wave  amplitude  at  a  receiver  location  r 
outside  of  S  is  given  by 


k  r 

u(r.<o)  =  — |dl  u(r’,<o)-^j-^(nr-+nr)«n, 


(2) 


where  n  is  the  unit  normal  of  the  closed  loop  S  (see  Figure  A1),  nr  is  the  unit 

vector  in  the  direction  of  propagation  of  the  surface  wave  incident  upon  S,  and 

r-r' 

nr  is  the  unit  vector  in  the  direction  of  the  receiver  from  the  loop,  nr  =  j^. 

The  following  derivation  of  Equation  (2)  for  the  surface  wave  ''2D" 
geometry  closely  follows  the  derivation  for  the  3D  scalar  wave  in  Klein  (1970). 
The  Fresnel-Kirchoff  integral  (Eq  2)  follows  for  any  field  that  obeys  the 
Helmholtz  wave  equation.  We  assume  that  two  fields  O  and  'V  each  obey  the 
Helmholtz  equation  over  a  surface  £, 
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I 


(3a) 

(3b) 


v2<i>-.k2<D  =  0. 
v2'P-k2'P  =  5(r). 


We  apply  Green's  theorem  to  the  two  fields, 

S  s  ^ 

where  the  surface  £  is  bounded  by  the  contour  S  formed  by  the  two  closed 
loops  Souter  and  Sinner  (see  Rgure  A2a).  The  area  integral  is  zero  since  ^  and 
O  are  solutions  to  Equation  (3). 


Sinner  Sourer 


We  choose  Smner  to  be  a  circle  of  radius  e  around  the  point  r  with  normal 
directed  toward  r.  r )  is  the  Green's  function  for  a  source  at  r,  and  the  integral 
over  Sinner  in  the  limit  that  e  goes  to  zero  becomes,  With  the  result  that 

for  k  |r-r'|  «  1 . 


Souter  Souter^ 


pi  (6) 


If  O  has  the  form  of  a  waveform  incident  on  the  loop  Souter  then  the  normal 
derivative  is  approximated  by  the  form  of 


dn 


-ikn*nr^(r'). 


(7) 


where  n^  is  the  direction  of  propagation  of  t&Cr').  Therefore,  we  have  the 
Fresnel-Kirchoff  integral  for  the  field  at  r  given  the  field  on  a  closed  loop  S 
enclosing  r, 
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Figure  A2.  The  region  £  is  bounded  by  Sinner  and  closed  loop  Souter  (A2a 
above).  Below,  the  loop  Souter  is  deformed  around  the  source 
and  scattering  region  with  a  section  So  that  goes  to  infinity  such 
that  Equation  (9)  becomes  Equation  (2)  with  a  reversal  in  the  sign 
of  the  normal  on  Souter- 
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(8) 


(nr*n-nrn). 


In  order  to  get  the  form  of  Equation  (2)  that  is  used  for  diffraction  calculations  we 
deform  the  outer  loop  to  infinity  and  surround  the  source  region  by  a  loop  as 
seen  in  Figure  A2b  connected  to  infinity  by  the  segment  Sq.  The  contribution 
from  the  So  segment  is  zero.  Note  that  this  is  equivalent  to  Figure  1  with  the 
exception  of  the  sign  of  the  normal  on  the  closed  loop  surrounding  the  source 
region  and  equation  (8)  becomes  equation  (2). 
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